Discontinuity of the chemical potential in RDMFT for open-shell systems 
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We employ reduced density-matrix lunctional theory in the calculation of the fundamental gap of 
open-shell systems. The formula for the calculation of the fundamental gap is derived with special 
attention to the spin of the neutral and the charged systems. We discuss the effects of different 
functionals as well as the changes due to different basis sets. Also, we investigate the importance of 
varying the natural orbitals for the calculation of the fundamental gap. 
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I. INTRODUCTION 

Density functional theory (DFT) 1 ' 2 is a powerful tool 
to calculate the electronic structure of atoms, molecules, 
and solids. Within DFT observables are given as func- 
tionals of the particle density. In reduced density-matrix 
functional theory (RDMFT) the 1-body reduced density 
matrix (1-RDM) is used as the basic variable. RDMFT 
is based on Gilbert's theorem 3 which proves that each 
ground-state observable can, in principle, be written as a 
functional of the 1-RDM. First-generation functionals 4-6 
perform very well in the description of the dissociation 
of small molecules. Second generation functionals were 
introduced very recently 7-9 which improved both the per- 
formance for small molecules 7-9 and also for the homo- 
geneous electron gas 10 . 

A key quantity in electronic structure calculations is 
the band gap for semiconductors and insulators. It is 
defined as the difference between the ionization potential 
/ and the electron affinity A 



that the fundamental gap is exactly given by 



where 



A = I — A, 



I = E tot (N - 1) - E tot (N) , 
A = E tot (N) - E tot (N + 1) . 



(1) 



(2) 
(3) 



Etot(N) denotes the ground-state energy of an N- 
electron system. In the chemistry literature A/2 is called 
the chemical hardness if the system is finite. For sim- 
plicity we use the term fundamental gap for both finite 
and extended systems throughout this article. We like 
to point out that the fundamental gap differs from what 
is known as the optical gap. The optical gap is given 
as the energy necessary to excite the system from the 
ground state to the first excited state. Therefore, its size 
is reduced by the binding energy of the created exciton 
compared to the fundamental gap. 

Within density functional theory it can be shown 11,12 



A = A 



KS 



A, 



(4) 



where Aks is the energy difference between the low- 
est unoccupied and the highest occupied Kohn-Sham 
states and A xc is the discontinuity of the exchange- 
correlation potential upon adding and subtracting a frac- 
tional charge. This discontinuity is zero for LDA and 
GGA, so Aks is the prediction for the gap within 
these approximations. However, this prediction deviates 
strongly from the experimental values. For semiconduc- 
tors the calculated gap underestimates the experimental 
value by typically 50%. In extreme cases, such as germa- 
nium, the gap vanishes within LDA. Interestingly, Aks 
for the exact-exchange functional is very close to the ex- 
perimental gap for several systems 13,14 . Unfortunately, 
in the case of exact exchange A xc is not zero and, in 
fact, was found to be much larger than Aks- Thus, if 
properly calculated, the band gaps within exact exchange 
are highly overestimated compared to the experimental 
values 13-16 . Exact exchange combined with RPA corre- 
lation was recently shown to yield results very close to 
the experimental values for Si, LiF, and solid Ar 16 (pro- 
vided the discontinuity A xc is properly included). Fi- 
nally, a recently introduced hybrid functional (HSE) 17,18 
is reported to give gaps in satisfactory agreement with 
experimental values for a set of 40 simple and binary 
semiconductors 19 . Especially, germanium is predicted a 
semiconductor with a gap of 0.56 eV. 

An alternative formula to (4) for the fundamental gap 
in DFT reads 20 



A= lim (u(N + ri) - (Jt(N - ri)) , 



(5) 



where \i is the chemical potential, and N is the particle 
number of the system. As Eq. (5) suggests, the chemical 
potential has a discontinuity at integer particle number 
N. In a recent paper 26 , we presented the analogous equa- 
tion within reduced density-matrix functional theory. In 
particular, we proved that the Lagrange multiplier used 



2 



to enforce the conservation of particle number is equal 
to the chemical potential. This theoretical development 
was applied to small finite and prototype periodic sys- 
tems with very promising results. We like to emphasize 
that the analogy between DFT and RDMFT is not at 
all trivial because of the ./V-representability condition in 
RDMFT. The occupation numbers are restricted to the 
interval [0, 1] which leads to border minima. For this rea- 
son the generalization of the proof of Eq. (5) from DFT 
to RDMFT is not straightforward. 

In the present work, we deduce a relationship similar to 
Eq. (5) for open-shell systems. The difficulty in generaliz- 
ing Eq. (5) to the open-shell case arizes from the fact that 
adding/subracting a spin-up electron to/from an open- 
shell ground state is not equivalent to adding/subtracting 
a spin-down electron. Open-shell systems were recently 
addressed in Ref. 21 where it was demonstrated that it is 
reasonable to introduce two Lagrange multipliers to keep 
the number of electrons in each spin channel fixed seper- 
ately. An alternative description of open-shell systems 
was introduced by Leiva and Piris 22 . In that desription, 
however, spin-up and spin-down occupations are equal 
for all orbitals except the open-shell ones which are fully 
occupied by the majority spin. The Lagrange multiplier 
is then spin independent. Here, we employ the treat- 
ment suggested in Ref. 21 where each of the two La- 
grange multipliers is a function of the two particle num- 
bers corresponding to the two spin components. In the 
present work, these particle numbers are assumed to be 
fractional. We show that a proper extension of Eq. (5) is 
possible with the resulting equation involving the discon- 
tinuities of both Lagrange multipliers. The derivation is 
presented in Section II. Section III contains results for a 
set of open-shell atoms and a comparison of the closed- 
and open-shell treatment for systems where the neutral 
system is actually closed-shell. We also investigate the 
performance of different functionals in the calculation of 
the fundamental gap. 



II. THE FUNDAMENTAL GAP IN RDMFT 

Rcduced-density-matrix-functional theory (RDMFT) 
uses the one-body reduced density matrix (1-RDM) 

7(x,x') = N J dx 2 ...dxjv**(x',X2, ...xjv)*(x,x 2 , ...xjy), 

(6) 

where \& denotes the many-body wave function and x = 
(r, a). Integration over dx means integration over space 
and summation over spin. Throughout this article we 
restrict ourselves, for simplicity, to the "collinear" case 
where 7(x, x') = 7(rer, r'er') is diagonal in spin space, i.e. 

7 (m,rV) = <W 7 CT (r,r'). (7) 



By diagonalizing 7 CT (r, r') one obtains the natural orbitals 
tp j<7 and the occupation numbers i.e. 

oo 

7 CT (r,r')=^n JCT ^ a (r')^(r). (8) 

.7=1 

To ensure the iV-representability of 7 the occupation 
numbers are restricted to the interval [0, 1] and sum up to 
the total number of particles N. In closed-shell systems 
the two spin-directions are identical, i.e. 

n rt = n 3U ( 9 ) 
¥>iT = fjl- (10) 

Within the spin-dependent formalism one can define 
spin-dependent electron affinities and ionization poten- 
tials by adding or removing an electron with specific spin 

r = E tot (N° -l,N*)-E tot (N°,N*), (11) 
A" = E tot {N° ,N S ) - E tot {N° + l,N s ) . (12) 

Here, E to t{N' T , N") representes the ground-state energy 
of a system with N = N a + N a electrons where N a is the 
number of electrons with spin a and N a is the number 
of electrons with the opposite spin, a. Consequently, 
the ionization potential and electron affinity defined in 
Eq. (3) are given by 

I = min{l\l 1 } , (13) 

A = m&x{A\ A 1 } , (14) 

cr 

i.e. they are respectively the smallest necessary energy 
for taking away an electron and the maximum energy 
gained by adding an electron to the neutral system. The 
fundamental gap then reads 

A = mm{f 1 I l }-m&x{A\A 1 }. (15) 

a a 

In order to derive a formula analogous to Eq. (5) for 
the fundamental gap (1) within RDMFT we follow the 
same path as in DFT 12 ' 23 ' 24 and extend the definition of 
the total-energy functional £aot[7] to systems with frac- 
tional particle number M. Throughout this paper we use 
the convention that N denotes an integer number of par- 
ticles and M a fractional. Such systems can be described 
as an ensemble consisting of an N- and an (iV+l)-particle 
state for N < M < N + 1. Let $n*,n* denote an TV- 
particle wave function with N = N a + N a where, as 
before, N" is the number of electrons with spin a and 
N 17 the number of particles with the opposite spin, a. 
We consider an ensemble where, compared to the charge- 
neutral (N a , N a ) system, the number of spin-er particles 
is increased by if*. The statistical operator describing 
such an ensemble is given by 

Dn°+7]° ,N s = (1 — V**) I ^N",N S ){^N",N' | 

+ lf |*JV- + 1,JV-)(*JV- + 1,JV- | • (16) 
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The expectation value of an operator O is then given by 
O = tr {Dn"+ V ",N'0^ . (17) 

In particular, for O = 7 CTl (r,r'), i.e. the operator repre- 
senting the 1-RDM of spin-di particles, we obtain 

7^ +7) „ >JV ,(r,r') = (1 - ?n7^(r,r') 

+ rfj N . +hN *(r,r'), (18) 

and for O = H, i.e. the Hamiltonian, we get the total 
ensemble energy 

E tot {N* + n a ,N*) = {l-r) a )E tot (N a ,N*) 

+ j ] ' J E tot (N a + l,N°). (19) 

We note in passing that the ensemble weights in Eq. (16) 
are such that the correct normalization of spin-up and 
spin-down densities is achieved, i.e. 



tfV 7 ^ +jrj7V5 (r, r) 
d 3 r^a +v „ N a (r, r) 
Reformulating (19) one obtains 



N° 
N° 



(20) 
(21) 



E tot {M°,N s )=E tot {N°,N s ) 

+ rf [E tot {N° + 1, N*) - E tot {N°, N 9 )] (22) 

for N a < M a = N a + rf < N a + 1. In analogy, for 
N a - 1 < M" = N a - 1 + rf < N a the total energy is 
given by 

E tot (M a , N°) = E tot {N° - 1, N s ) 

+ rf [E tot (N°, N°) - E tot {N° - 1, N s )} . (23) 

In other words, the total energy depends linearly on 77 " 
with slope -A a for N a < M" <N a + l and slope -T 
for N" - 1 < M a < N a . Since A" and I" are in general 
not the same, the derivative dE tot (M er ,N^)/dM (T has a 
discontinuity at integer particle number N a . From Eqs. 
(11)-(15), one can conclude that the fundamental gap is 
given by 



A = min < lim 



8E tot {M\Mt) 



max < lim 

a J7 "— >0+ 



dE tot {M\M^) 



N"+rt" ,N° 



dM c 



(24) 



In Ref 21, we argued that, for open-shell systems, the 
following functional should be minimized 



F[j}=E tot [y] 



I E n H ~ M? 



(25) 



The Lagrange multipliers fi' and fi* are introduced to 
achieve given particle numbers and MK To prove 
the formula for the fundamental gap we first show that 
these Lagrange multipliers are nothing but the chemical 
potentials, i.e. 



»° { MiM) = dEtot{ f ,Mi) 



mJ ,m, 1 



(26) 



The derivation of this formula differs significantly from 
the derivation of its counterpart in DFT due to the above 
mentioned ./V-representability constraint. In order for 
the 1-RDM to be connected to an anti-symmetric TV- 
particle wave function its occupation numbers have to 
be restricted to the interval [0,1] 25 . One can show that 
the same constraint ensures ensemble A^-representability 
for fractional particle number. As a result of this addi- 
tional constraint, 5F/Sj need not vanish at the minimum 
energy. It is possible that certain occupation numbers 
are pinned at the border of the interval while the true 
minimum is obtained for values of rij C outside this in- 
terval. The functional F then has a border minimum, 
and therefore non-vanishing derivative, in all directions 
where occupation numbers arc pinned at zero or one. 
We investigate the difference 

E tot (M a + rf, M°) - E tot (M°, M°) = 

E YlM'+rf,M a ] ~ E [lM",M*] • (27) 

A Taylor expansion of E["/M<*+r]<' ,m s ] around "/m°,m s 
yields 

E tot (M° + rf, M°) - E tot (M°, M°) = 
£ ffd 3 rd 3 r' 6Etot 



£7 CTl (r,r') ^ 
x (7^ + ^, M 4r,r')-7^,M^(r,r')). (28) 
For the functional derivative we employ (25) and obtain 
SE tot SF 



S^o-i ( r ,r') (57 CTl (r,r') 

00 e 

5n n 



T E 



/-E 



5n 



j I 



J (57 CT i(r,r') j-^ d^" 1 ( r > r 



(29) 



The first term on the right is evaluated via the functional 
chain rule, i.e. 



SE t , 



d"f ai (r, r') 



n = E E 



Tt i=i 



/ 6F 



/ 



d 3 r" 



SF 6<p 
5<p jcr (T") 8Y 1 (r,r') 



+ c.c 



(30) 



J =1 



At the solution point, the variation with respect to the 
natural orbitals vanishes such that the second term on 
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the right is zero. The variation with respect to the oc- 
cupation numbers, however, need not vanish due to the 
7V-representability constraint. Equation (28) therefore 
reduces to 

E to t{M° + rf, M°) - E tot {M°, M s ) = 



EE 

<T1 = U P 



d\d^l— ^ 



8n pai <S7°"i(r,r') 

X [7A/"+r;",M s ( r ' r ) — 7A/",A/ S ( r > r ) 

EE// d"rd\'^ 



faff! f r r A 

where the first sum runs only over those occupation num- 
bers pinned to the border of the interval. The variation of 



the occupation numbers can be calculated applying first 
order perturbation theory to the eigenvalue equation of 
the 1-RDM 



d s rV r (r>r>*r( r/ ) = nj*<Pi*(*) ( 32 ) 



which yields 



8^ a (v, r') 



(33) 



In addition, we write the eigenvalues and eigenfunctions 



of 



"l 'M" +7J" ,A/ d 



as 



M°+rf ,M" . r M"+rf ,W 



(34) 

where (pj ai and rij ai denote the natural orbitals and oc- 
cupation numbers of 7^ M &. Equation (31) then re- 
duces to 



J 




where we only kept terms up to first order. In order for 
the natural orbitals to remain normalized the changes 
have to be orthogonal to the original orbitals, i.e. 



i ry* rT (r)8y jrT {r)=Q. 



(36) 



Therefore, the integrals one the right-hand-side of Eq. 
(35) vanish. The sum over all changes in the occupation 
numbers has to give rf in order for the new occupa- 
tion numbers to sum up to the correct particle number. 
Hence, we obtain 

E tot (M° + rf, M°) - E tot (M a , M°) = 

^ + E E ( 37 ) 



<Tl=U P 



Finally, we discuss the contribution of the pinned states. 
As stated before, for these states 5F/5n p is different from 
zero and the true minimum of the functional lies outside 
the interval [0,1]. More specifically, it lies at a finite 
distance from the border of the interval such that the 
addition or subtraction of an infinitesimal fraction rf of 
a particle cannot move the minimum into the interval. 
Therefore, these particle numbers remain pinned upon 
adding or subtracting an infinitesimal rf , i.e. 8n pai is 



zero in the limit rf — > 0. We therefore conclude 
if {Ml, = 



lim 

77°"— >0+ 



Etot(M° + rf, M g ) - E tot {M°, M s ) 
dE tot (M\M-t) 



Ml, Ml 



4, Ml 



(38) 



Using Equation (24) we obtain the final result for the 
fundamental gap 



A = min [ lim fi a (M a + rf , M a ) 
-max( lim /i"(M° 



rf t M')\. (39) 



The derivation of Eq. (39) concerns the exact exchange- 
correlation energy functional of the 1-RDM. Since only 
approximations are available, the question is whether 
Eq. (39) is still useful. This question is the main sub- 
ject of the next section. 

In Ref. 26, a single, spin independent fi (for closed- 
shell systems) was shown to have a discontinuity as a 
function of a fractional total number of electrons which 
is equally distributed in the two spin channels. The ap- 
plication of that theory to an open-shell system would 
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give the spin resolved \i" as a function of a unique M. In 
the present work, we add/subtract a fractional part of an 
electron to/from a specific spin channel. Consequently, 
the system becomes an open-shell system even if the neu- 
tral system is closed-shell. Thus, we have four functions 
^(Ml,N l ), fi l (N\M l ), ^(N\M l ), and /^(M 1 ", N l ), 
where , arc fixed to the integer values of the neu- 
tral system. Of these four, only the first two show a 
discontinuity. The correct gap is then given by Eq. (39), 
where the min and max functions take care of the selec- 
tion of the smallest ionization potential and the largest 
electron affinity. Alternatively, one can employ Eqs. (1)- 
(3) for the calculation of the fundamental gap. Both 
approaches are exact, in the sense that, given the exact 
functional of 7, they both reproduce the fundamental 
gap. It is interesting to see if they give the same num- 
bers for approximate functionals as well. This is also one 
of the questions we address in the next section. 

To answer the above questions, one needs to minimize 
the approximate functionals for fractional number of par- 
ticles to get ^}'' L (M a , N a ). The extension of the mini- 
mization procedure to fractional particle numbers, which 
is in complete accordance with the proof we presented 
above, requires us to perform the minimization in the 
domain of 7^ r 1 , T+r;<T N s, which are given by Eq. (18). In 
principle one then has to minimize the total energy with 
respect to 7^ NS and 7^+1 at* under the known N- 
represcntability constraints that their occupation num- 
bers are between and 1 and sum up to the correct par- 
ticle numbers. However, this procedure, involving the 
density matrices for N and N + 1 particles, is not very 
practical. On the contrary, it is desirable to minimize 
with respect to 7^ r 1 „ + ^„ N a directly under the appropriate 
constraints. We prove elsewhere 27 that the appropriate 
constraints for such a minimization are 

0<ngf^<l Vi, £ngf^=M-. (40) 

3 

In other words, the domain of 7^ 1 CT+1)tT N „ which can be 
represented as the weighted average Eq. (18) is identical 
to the domain of "f^^v N & whose eigenvalues satisfy 
Eq. (40). The above statement is quite significant since 
the constraint of Eq. (40) is much simpler and completely 
analogous to the case of integer particle numbers. The 
implementation is therefore a rather simple extension of 
the case of integer particle numbers. 



III. NUMERICAL RESULTS 

In this section, we study the behavior of fj, as a func- 
tion of the fractional particle number for some atoms and 
molecules using approximate functionals of the 1-RDM. 
Our aim is to investigate whether there exists a discon- 
tinuity in fi(M) and how it compares to the fundamen- 
tal gap. The implementation we used for finite systems 
can be applied to both closed- and open-shell 21 configu- 



rations. Some results for closed-shell systems were pre- 
sented in Ref. 26. Here, we give an extended analysis for 
both closed- and open-shell systems. 

For the open-shell treatment, we use the extension 
of the functional of Goedecker-Umrigar 5 described in 
Ref. 21. We also investigate whether other functionals 
reproduce a discontinuity in a closed-shell treatment. For 
this purpose we consider the functionals of Piris 8 ' 9 , where 
the self-interaction (SI) terms are explicitly removed, and 
the Miillcr functional and the most recent BBC func- 
tionals of Gritscnko et al 7 which contain self-interaction 
terms. 

The implementation is based on the GAMESS 
program 28 which we use for the calculation of the one and 
two-electron integrals. The minimization with respect to 
the occupation numbers and natural orbitals is then per- 
formed using the conjugate gradient method. Our pro- 
gram treats both closed- as well as open-shell systems 
using the restricted open-shell RDMFT 21 . In short, we 
assume spin-dependent occupation numbers (and chemi- 
cal potentials) but spin-independent natural orbitals. In 
that way, our method is in complete analogy to spin re- 
stricted open-shell Hartree-Fock. 

In Fig. la, we show n(M) for the LiH molecule us- 
ing the GU functional in the closed-shell treatment, i.e. 
the extra charge is equally distributed over the two spin 
channels. Fig. lb shows ^ T (M T ,iV x ) and ^ X (M T ,A^) 
for the open-shell treatment of the LiH molecule, using 
again the GU functional. In the open-shell treatment the 
additional charge is exclusively added to one spin chan- 
nel, and here we choose the spin-up channel. Clearly, 
/i(M) in Fig. la and ^ {M\N^) in Fig. lb show a pro- 
nounced step which resembles the discontinuity that one 
expects for the exact functional. This step has two im- 
portant features: the first is that it occurs not exactly 
at M = 4, i.e. the exact, integer number of electrons. 
It is rather shifted slightly to the right. The shift is of 
the order of 0.05 of an electron in Fig. la and is reduced 
to 0.02 in Fig. lb. A closer look at the solution reveals 
that the bottom of the step appears exactly at the point 
where the occupation number of the HOMO gets equal 
to one. After that point it has to remain one due to the 
A^-represcntability constraints, Eq. (40). The pinning of 
the occupation number of the HOMO to one results in 
the rapid increase of /i. Since adding charge to one spin 
channel only results in faster pinning of the HOMO state 
it is not surprising that the step in the open-shell treat- 
ment is shifted to the left. Upon increasing the extra 
charge further, [i is a smooth function, i.e. the upper 
edge of the step is rounded off. In the closed-shell treat- 
ment fi(M) shows a linear dependence outside the step 
region, which is significantly reduced in fj,' in closer re- 
semblance to the exact behavior. A more detailed in- 
vestigation reveals that the slope of n(M) is the average 
of the slopes of ^(M\N l ) and n l (M\N l ). To ex- 
tract a value for the discontinuity we use a backwards 
projection as shown in Fig. 1. This method reduces to 
the exact discontinuity if [i is a true step function. The 
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FIG. 1: The behavior of fi as a function of a fractional electron 
number M for the LiH molecule in a closed-shell treatment (a) 
and fj?'+ (M T , iV^ ) for an open-shell treatment (b). For com- 
parison, the experimental and CI values of the fundamental 
gap are included. 



extracted values, as well as the gaps of other finite sys- 
tems are given in Table I. We should also keep in mind 
that DFT methods like LDA and GGA underestimate 
the gap by typically 50%. Although the procedure of 
backwards projection might seem rather crude and ar- 
bitrary, wc should mention that the agreement with ex- 
periment is rather satisfactory for both close- and open- 
shell treatment. As one can see, for LiH, the quantitative 
agreement is slightly better for a closed-shell treatment. 
Nevertheless, the open-shell treatment should be prefered 
because (i then resembles the exact step function much 
closer making the backward projection less ambigious. 

For open-shell systems, varying or M* is not equiv- 
alent anymore. Thus, we can study the behavior of both 
[V and fi* as functions of M' or MK We investigate 
the open-shell atoms Li, Na, and F varying M' or M* 
away from the neutral configurations. In the following, 
we use the convention that spin up is always the ma- 
jority spin channel. In Fig. 2, we show the results for 
fi a for the Li atom. Only the chemical potential cor- 
responding to the spin direction whose particle number 
is changed shows a discontinuity as already observed for 
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FIG. 2: The behavior of /i CT as a function of an electron 
fraction r]' 7 added (subtracted) to the neutral system for the 
Li atom. In the inset, we show an enlargement of the re- 
gion where we extract the value for the gap from the differ- 
ence of the upper level of ^(N^, M+) and the lower level of 
y}{M\N l ). 
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of the gap are extracted from. 
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System 


RDMFT 


RDMFT 


Other 


Experiment 




H(M) step Eqs. (l)-(3) theoretical 




Li 


0.18 


0.202 


0.175" 


0.175 6 


Na 


0.18 


0.198 


0.169 c 


0.169 6 


F 


0.54 


0.549 




0.514 6 


LiH 


0.27 d ,0.29 e 


0.271 


0.286 / 


0.271 9 



TABLE I: The prediction for the fundamental gap for sev- 
eral atoms and small molecules using the size of the step of 
/i(M), and a direct calculation through Eqs. (l)-(3) for the 
GU functional compared with experimental and other theo- 
retical values. For the direct application of the Eqs. (l)-(3), 
the total energies of the positive and negative ions were cal- 
culated. 

a QCI from Ref. 29 
b from Ref. 30 

c Ionization potential from 29 , electron affinity from 31 
d Closed-shell treatment 
E Open-shell treatment 

^ CISD using the same basis set as in RDMFT 

9 Ionization potential from 32 , electron affinity from 33 



the LiH molecule. Therefore, we only plot fj,'(M',N^) 
and /j,^(N^, M^). Again, pronounced steps resembling 
the discontinuity of the exact theory are present. The 
prediction for the gap is then selected using Eq. (39) and 
the backwards extrapolation procedure described earlier. 
The values obtained for the gaps are listed in Table I. 
According to Eq. (39), the gap for the Li atom is given 
by the difference between the backwards projected upper 
part of fi l (N\M l ) and the lower part of ^ {M\ N l ). 
In Fig. 3, we show the analogous results for the Na and 
F atoms. The picture for the Na atom is very similar to 
Li. On the other hand, for the F atom, the gap is given 
by fi^(N^ , M^) alone. It is interesting that the position 
of the upper and lower parts of the corresponds to 
the actual process of adding and removing electrons to 
the system. Thus, for Li and Na atoms, it is favorable to 
remove an electron from the majority spin channel (up) 
and add an extra electron to the minority spin channel 
(down). As a consequence the gap is given by the dif- 
ference between the upper part of n^(N^, M^-) and the 
lower part of /jJ(M', N^). For a F atom, on the other 
hand, it is favorable to add an electron to, or remove 
from, the minority spin channel. Thus, the gap is given 
by li 1 (N^,M 1 ) alone. 

In Table I we give the results obtained by the backward 
extrapolation for the systems discussed in this paper. As 
one can see, they agree very well with experimental val- 
ues for the fundamental gap as well as other theoretical 
calculations. For finite systems, one can also calculate 
the gap by performing three total energy calculations, 
for the N, the N + 1 and N — 1 particle systems and use 
Eqs. (1-3). The values for the gap obtained in this way 
are given in Table I for comparison. One should keep 
in mind that for solid state systems, this procedure does 
not apply because the addition or the removal of a sin- 
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FIG. 4: The function fj,(M) for the He atom using the cc- 
PVQZ basis set without and with an additional very diffuse 
s-type basis function. 



gle electron to an infinite solid is meaningless. For such 
systems, the recipe introduced in this work is expected 
to be valuable. 

Of course the question arises whether the system with 
excess charge is correctly described by the basis set we 
used. Usually, atomic basis sets are optimized to cor- 
rectly describe the neutral system resulting in basis func- 
tions which arc all localized. Therefore, the charged sys- 
tem might be predicted to have a localized bound state 
despite the fact that the configuration of a neutral atom 
and a free completely delocalized electron is energetically 
favorable. A prominent example of a system not having 
a negative ion is the He-atom. We study the behavior of 
/i(M) with two different basis-sets: the CC-PVQZ basis- 
set and CC-PVQZ enlarged by a very diffuse s-type func- 
tion. As one can see in Fig. 4, the state of the additional 
fractional electron is better described by the enlarged ba- 
sis set. In this case, the electron affinity is zero and the 
gap is given by the IP alone. Interestingly, the inclusion 
of a diffuse function leads to a sharper step of /x(M) in 
close resemblance to the discontinuity of the exact func- 
tional. We also add extra diffuse functions in the basis- 
sets of both Li and H in the calculation of /x(M) for the 
LiH molecule. We do not observe any effect on /u(M), 
which is a clear evidence for the fact that LiH binds an 
extra electron and that the localized basis-set is appro- 
priate for describing the state of the charged system. 

In order to investigate the importance of the variation 
of the natural orbitals for the discontinuity of /i we per- 
form, apart from the full variation described so far, a 
calculation where only the occupation numbers are opti- 
mized while for the natural orbitals we keep the initial 
Hartree-Fock orbitals. In Fig. 5 we compare these two 
procedures for both a closed- and an open-shell calcula- 
tion. As one can see from the plots, the main contribu- 
tion to the discontinuity arises from the variation of the 
occupation numbers. In the closed-shell calculation we 
obtain a discontinuity of 0.27 Ha for the full variation 
compared to 0.26 Ha if we vary the occupation numbers 
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FIG. 5: The function (i{M) (n 1 (M r , N 1 )) for the LiH 
molecule using the closed-shell (a) and the open-shell treat- 
ment (b), with occupation number variation (using the 
Hartree Fock orbitals) and with full variation (both occupa- 
tion numbers and the orbitals). 
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FIG. 6: The function n(M) for the LiH molecule using the 
closed-shell treatment for the Goedecker-Umrigar, the Miiller, 
BBC1, BBC2, BBC3, and PNOF functionals. The first and 
the last involve a complete removal of the SI terms. Only 
these two reproduce a pronounced step in resemblance to the 
discontinuity of the exact theory. 



Therefore, we conclude that the complete removal of the 
SI terms is essential for obtaining the correct behavior 
of fi(M). The size of the step of fi(M) for the PNOF 
is 0.30 Ha and compares well with experiment (see Ta- 
ble I) . As a test, we also tried a modified version of BBC3 
where we removed the SI terms completely. Consistent 
with the above conclusion, it also produces a step which 
is almost identical to PNOF. Additionally, this modified 
BBC3, like the GU functional, gives an accurate measure 
of the correlation energy at the equilibrium distance, but 
fails completely at the dissociation limit. 



only. In other words, only about 4% of the discontinu- 
ity are due to the optimization of the natural orbitals. 
This picture remains unchanged if we use the open-shell 
procedure where we obtain 0.31 Ha for the full variation 
and 0.29 Ha for the variation of the occupation numbers 
alone. 

In all the calculations presented so far, we have used 
the functional of Goedecker and Umrigar, which involves 
the complete removal of the self-interaction terms. It is 
interesting to study the behavior of /u(M) using different 
functionals, like for instance the recent BBC function- 
als of Gritscnko et al 7 and the PNOF of Piris 8 ' 9 . In the 
BBC1 and BBC2 functionals, the SI terms are present 
while in the BBC3, they are partially removed. How- 
ever, the SI terms for the bonding and the anti-bonding 
orbitals remain. In the PNOF they are fully removed as 
in GU. In Fig. 6, we plot fi(M) for LiH using the closed- 
shell treatment, for all these functionals. Surprisingly, 
only GU and PNOF show a pronounced step which com- 
pares well with the fundamental gap. The other func- 
tionals show either a completely smooth behavior or, in 
the case of BBC3, a small kink in the wrong direction. 



IV. CONCLUSION 

We have presented a formalism to calculate the fun- 
damental gap within RDMFT for both open- as well as 
closed-shell systems. Our numerical results show that 
even for systems where the neutral system is closed-shell 
the results for the chemical potential are closer to the ex- 
act step function if an open-shell treatment is employed 
because adding charge of a specific spin to the system 
makes it open-shell. The application to several open-shell 
systems gives a very good agreement with experimental 
values in all cases. Also, the steps in the chemical po- 
tentials are such that they resemble the spin dependence 
of the ionization potential and the electron affinity of 
the real system. Our investigation of a possible basis set 
dependence reveals that it is necessary to include very 
diffuse states in the basis set in case the system does not 
bind extra charge. Whenever the system does bind ex- 
tra charge the results are independent of the inclusion 
of the diffuse state in the basis set. To estimate the 
contribution of the occupation numbers and the natural 
orbitals to the fundamental gap we compared the results 
for the LiH molecule using a full variation and a varia- 
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tion of the occupation numbers only. We found that over 
90% of the fundamental gap are due to the occupation 
numbers. This finding was confirmed for several other 
systems so far and we believe that it shows a general 
feature of RDMFT calculations. Finally, we investigated 
the behavior of several different functionals for the cal- 
culation of the fundamental gap. From our results we 
conclude that the exclusion of the self-interaction for all 
natural orbitals is essential to obtain reasonable results. 
Functionals without any removal of self-interaction sim- 
ply yield a continuous chemical potential. 

The present work is a contribution to the subject 
of calculating the fundamental gap of materials within 
RDMFT. The hope is that this theory gives results closer 
to experiment than DFT for this fundamental problem. 



It is our belief that the theoretical development presented 
in this work will have a significant impact in the appli- 
cation of RDMFT to periodic systems. 
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